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1. INTRODUCTION 

This paper describes some developments in fast elliptic solu- 
tion algorithms and methods for their application in computa- 
tional fluid dynamics . The description concerns mostly past 
work, but includes brief remarks on work in progress and some 
expectations for the future of similar methods. 

The primary importance of the fast, direct elliptic 
solvers lies in the possibilities for their application in 
complex practical problems in fluid dynamics as well as in 
other areas of computational physics. The work summarized 
here is motivated by the need for the efficient solution of 
practical problems. 

The partial differential equations governing fluid- 
dynamic problems are generally nonlinear. The approach empha- 
sized here for using direct elliptic solvers in the numerical 
solution of nonlinear problems is referred to as a seraidirect 
method, which is a globally implicit iteration scheme that 
uses a direct elliptic solver as the driving algorithm of the 
iteration. 

The following topics are discussed in connection with the 
semidirect method: (a) the equations of fluid dynamics that 

have been treated or that are expected to be treated, (b) the 
use of direct solvers for the globally implicit iteration, 

(c) a Cauchy-Riemann solver and its relation to Poisson 


* 

Contribution to the GAMM Workshop on Fast Solution 
Methods for the Discretized Poisson Equation, held at the 
Kernf orschungszentrum, Karlsruhe, Germany, March 3-4, 1977. 
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problems, (d) some developments needed for application of the 
semidirect method to the mixed elliptic-hyperbolic equations 
of steady, inviscid transonic flow, and (e) the treatment of 
interior boundary conditions. 

2. SEMIDIRECT ITERATIVE METHODS 

Consider the equations representing four levels of approxima- 
tion of fluid flows: (I) Navier-Stokes equations, (II) Euler 

equations, (III) Potential-flow equations, and (IV) Small- 
disturbance-potential-flow equations. The Navier-Stokes 
equations (I) describe a large class of fluid flows. However, 
when the transport processes including viscosity and heat 
conduction are negligible, the Euler equations (II) aovern the 
"inviscid flow" adequately. 1 -f there are no sources of 
rotationality , such as strong curved shock waves, in the 
inviscid flow, the potential-flow approximation (III) is 
valid. The equations are further simplified (IV) if only 
small disturbances affect the flow. Fromm /35/ was the first 
to apply a direct elliptic solver in nonlinear fluid dynamics. 
He used Buneman's /12/ Poisson solver to treat part of the 
Navier-Stokes equations for an unsteady, nearly incompressible, 
laminar flow. Since then, many others have used Poisson 
solvers in similar ways. Roache /97/ then used a Poisson 
solver in a different way, as the total driving algorithm in 
an iterative scheme (i.e., a semidirect method), to solve 
the Navier-Stokes equations for a simple viscous flow. At 
about the same time, Widlund /141/ and Concus and Golub /18/ 
were exploring similar iterative methods for nonseparable 
(but linear) elliptic equations. 

To further describe the semidirect method, let us con- 
sider the equations of steady potential flow (III) , 

U x + V y = - (1/P) <UP X + V Py ) , U y - V x = 0 , ( 2 . la ,b) 

p = F(p 2 U 2 +p 2 V 2 ) = [ 1+ (1/2) (T-l)M^(l-U 2 -V 2 ) ] !/ (T-1 ) (2.1c) 

where p is the dimensionless mass density, U and V are the 
dimensionless Cartesian velocity components, and y and M m 
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are constants. An equivalent system for the stream function 
4' defined by pU = ^ , pv = is 

* xx + V = (1/p) ^A + PyV ■ p - >?(* x + *y) • < 2 ' 2a ' b > 

Another equivalent system can be obtained from eqs . (2.1) 

for a velocity potential $ defined by U = 4> x , V = ‘fcy 

If the right side of either eq. (2.1a) or (2.2a) were 
known, these would be just either nonhomogeneous Cauchy- 
Riemann equations or a Poisson equation. A simple form of 
the semidirect iterative method treats the right side as 
known from a previous iteration and uses a direct elliptic 
solver for the left side; the iteration is continued until 
the solution converges. This method has been used (in a 
preliminary version of /73/*) to solve a finite-difference 
form of eqs. (2.2) for purely subsonic flow over lifting air- 
foils. In a similar way, one could use a direct Cauchy- 
Riemann solver /61, 74,71/ to solve instead eqs. (2.1). In 
these subsonic applications, central differences are used 
everywhere because the equations are elliptic. 

In steady, inviscid transonic flow, on the other hand, 
the equations are elliptic in the subsonic portions of the 
flow, but hyperbolic in the supersonic portions. In addi- 
tion, discontinuities in the variables must be allowed for 
because the transition from a supersonic region to a subsonic 
region may occur through a shock wave. Jameson /53/ has 
solved the so-called 11 full-potential equation" (III) in terms 
of the velocity potential for transonic flows using a line- 
relaxation iterative method (with a direct Poisson solver 
used to accelerate convergence, but not as the driving 
algorithm). Work is currently in progress to solve eqs. (2.1) 
for steady transonic flow over lifting airfoils using a form 
of the semidirect method. 

The semidirect method has been used /73, 67-69/ to Solve 
the transonic-small-disturbance equations (IV) 4 


*AIAA Paper 74-11, Jan. 1974. 
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U x + V y = M « u x + 7 ^ +1 ) m L uu x * u y " v x = 0 * (2 . 3a ,b) 

where the perturbation velocities u and v are given by 

U = 1 + tu , V = rv , (2.4) 

r is a small parameter (e.g. , airfoil thickness ratio) , 
and M m is the free-stream Mach number. Eqs . (2.3) are 

locally elliptic or hyperbolic depending on whether u > u*, 
the critical value. If the flow is "subcritical" the semidi- 
rect method with central differences is straightforward. In 
supercritical cases, the hyperbolic regions and shock -wave 
discontinuities require special treatment; the additional 
developments needed are described in section 4 . 

3. DIRECT CAUCHY- RIEMANN SOLVERS 

A direct solution algorithm for the elliptic equations 

u + v = s(x,y) , u - v = -w(x,y) , (3.1a,b) 

x y A 

is called a Cauchy-Riemann solver ,'61,74,71/. (The earliest 
"direct" Cauchy-Riemann solver appears to have been given by 
Gates and von Rosenberg /3 7/', but the more recent "fast" 
methods do not have some of the apparent limitations of the 
approach used in /37/. ) For the applications, additional 
terms can be added to the left side of eqs . (3.1) /74,71/. 

If = 0 everywhere, one can define a function 0 

(e.g., velocity potential) by u = 0 , v = 0 , and obtain a 

x y 

Poisson equation for 0. Or if s = 0 everywhere, one can 
define a function 0 (e.g., stream function) by u = 0y r 

v = -0 , and obtain a Poisson equation for 0. The 

A 

0-formulation allows point sources (s) , but not point vor- 
tices. The 0-formulation allows point vortices (cj) , but not 
point sources. But the u-v formulation allows both s and 
to to be nonzero. Furthermore, conditions such as at solid 
boundaries or across shock waves in the semidirect method 
&XO 8\M)V v i.fied most naturally in terms of u and v. For a 
given problem, the boundary conditions may make the Cauchy- 
Rit »\&hn formulation preferable over a Poisson formulation. 
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For given Poisson problems, one can solve equivalent 
Cauchy-Riemann problems and, in the process, obtain the 
Poisson solution. To illustrate this while formulating the 
Cauchy-Riemann algorithm /71/, let us consider the discret- 
ized problems as follows. Let the indices for u be i and 

j ' , the indices for v be 


r+i- 
j - 
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j-i 
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u 
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Fig. 1.- Indices for discret- 
ization on staggered grid. 


i' and j, and let the points 
where u and v are defined 
be oriented in a staggered 
grid as in Fig. 1. The 
"continuity equation" (3.1a), 
central differenced about the 
point and the "rota- 

tionality equation" (3.1b), 
central differenced about 
i,jr are 


(u i,j' - u i-l , j ' )/Ax + <V,j - V i',;j-1 ,/Ay “ . (3-2a) 


- U i,j' )/Ay - - V i',j)/ Ax = -“i,j • < 3 - 2b > 


On the computation field in Fig, 2, u is defined at points on 
the left and right boundaries and v is defined on the top 
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Fig. 2.- Meshes for 
Poisson-Dirichlet , 
Poisson-Neumann , and 
equivalent Cauchy- 
Riemann problems. 


E. D. Martin 


6 


and bottom boundaries. With this mesh, boundary values of u 
and v are specified respectively on those boundary segments. 
Equation (3.2a) is to be satisfied at all interior points 
designated by <f> on Fig. 2, and eq. (3.2b) is to be satis- 
fied at all interior points designated by Let us relate 

this configuration and Cauchy-Riemann (CR) problem to both 
the Poisson-Dirichlet (P-D) , or first boundary-value problem 
(n^p = 1) , and the Poisson-Neumann (P-N) , or second boundary- 
value problem (n bvp = 2 ). if all s ± , , - 0 , define 


u. . , = (i p. . 
1,3 i r 3 




V . . . = ( ^ . 

i',3 lv i-i,3 


- 0. .)/Ax 

X / J 


(3. 3a,b) 

Equations (3.2) are then equivalent to a discretized Poisson 

equation for , with values of & specified on the boundary 

(n bvp = 1). Thus, a P-D problem can be solved by using an 

algorithm for eqs. (3.2) on the mesh of Fig. 2, and the 

values of $ . . are then obtained from eqs. (3.3). If, on 

i , 3 

the other hand, all oi. - o, define 

1 r 3 


U, 


1 r 3 


V . , . 

1 1 ,3 


^i' + i, j ' " *i',j ,)/Mx ‘ 
= <*i' f j'+i “ *i',j' )/Ay 


(3 . 4a,b) 


Eqs. (3.2) are then equivalent to a discretized Poisson equa- 
tion for 9 , with normal differences of $ centered on each 
boundary; and specifying u and v on the boundaries is equiv- 
alent to specifying the normal differences of <j> (n bvp = 2 ) . 

A P-N problem can therefore be solved by using the algorithm 
for eqs. (3.2) on the mesh of Fig. 2 (with the boundary 
values of u and v specified in the same configuration as 

for n, = 1), and 9., is obtained from eqs. (3.4). 

bvp x ,3 

With an algorithm for solving eqs. (3.2) on Fig. 2, 

there is no essential difference in the treatments for 

n, =1 and 2 except that in the one case all s . , . , 

bvp x ,3 

zero and (3.3) gives the final solution, whereas in the 
other case all co. . are zero and (3.4) gives the final 
solution. In previous formulations of Poisson solvers, a 


are 
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different coefficient matrix results for the two problems, 
and the P-N solution has usually been more time consuming. 

When eqs. (3.2) are written for all i,j, an equation 
is obtained /71/ that reduces to the very regular block- 
tridiagonal matrix equation 

TV = F , (3.5) 

in which T is a block- tridiagonal matrix having entries 

(-1, C, -I) , £ is a unit matrix, C is a scalar-tridiagonal 

matrix having entries (-o', @ -a) , a = (Ay/Ax) 2 , P ^ = 2 + a 

for the first and last diagonal elements, and otherwise 

P ^ = 2 + 2a . The column vector V is composed of all the 

v., . and the column vector F has elements that are simple 

i 1 3 

combinations of the s . , . , and to. 

i'/3’ t/3 

Note that (a) all the u. have been eliminated from 

1/3 

the equation to be solved, (b) eq . (3.5) can be solved by 

any fast block-tridiagonal-equation solver, (c) after all 
the v. , . are known, one can obtain the u. . , quickly 

^* / J , r D 

from eq. (3.2a) for j' = n T and from eq. (3.2b) for all 
other j'. If only a Poisson solution is wanted, one need 
not determine all u^ . , but may simply find either the P-D 
or the P-N solution from eq. (3.3b) or (3.4b). 


A FORTRAN computer program (UVSOLV) was written for 
testing this algorithm and for comparison with other fast 
Poisson solvers. The "standard problem" for the comparison 
was the discretized Poisson equation with either Dirichlet 
^ n bvp = ^ or Neumann ( n bvp “ conditions. That problem 
was equivalent to eqs. (3.2) with either (3.3) or (3.4) and 


S i 1 , j ' ^ n bvp -1 ^ ^k ,£^i ' j ' ' 


to. . = (n, -2) A. 0 \p. . 

i,D bvp k,«*i,3 


A k = 2[ cos (ka) + cos (2b) -2} , a = n/m 1 , b = v/n 1 


(3. 6a, b) 
(3.7) 


where the functions put into eqs. (3.6) are the exact 
discretized solutions for n b V p = 1 anc ^ 2 / respectively, 


On the vertical boundaries, i = 0 and ; on uhe 
horizontal, j = 0 and n^ 1 
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■ - sin(ika) sin(jfb) , (3.8a) 

1 / J 

0., . , = cost (i 1 - l/2)ka] cost (j 1 - 1/ 2 ) fib] . (3.8b) 

1 t 3 

For both cases, the integer index i (or i') ranges from 0 

to M + 1 and j (or j') ranges from 0 to N + 1 where 

M = m + n. - 2, N = n + n, -2. Zero boundary values 

i bvp i bvp J 

for u and v are specified. 

The computer program solved eq . (3.5) by recursive 

cyclic reduction, with the inner loops (scalar tridiagonal 
solutions) done by Gaussian elimination. The program was 
run on the Control Data 7600 computer at Ames Research Cen- 
ter, and a slightly modi.fied program was then run on the IBM 
370/168 at the Kernf orschunaszentrum at Karlsruhe, Germany. 
The test cases used nij = 128, n ( = 32 for both the P-D and 
P-N problems. For the 128X32 mesh, the CPU time on the 7600 
(FTN 4.5 compiler) for both the P-D and P-N problems was 
0.054 (±0.001) sec. On the 370/168 (H compiler), the solu- 
tion times were 0.18 (xO.Ol) sec and 0.19 (±0.01) sec. 

Different values of k and £ were used, some of which 
provided a severe test of the accuracy. Single-precision 
truncation errors on the IBM machine were excessive for a few 
cases. Of 18 cases, 6 had maximum errors larger than 10 -3 , 
one larger than 10 -2 (6.64X10 -2 ). (A subsequent modification 
reduced the largest error to 3X10~ 3 . ) However, on the 7600, 
all errors were less than 10 -9 , with all but two being less 
than 10 -1 0 . On an IBM 360/67 at Ames Research Center (having 
the same precision as the 370/168) no accuracy difficulties 
have been encountered in applications to physical problems. 

It has also been pointed out that a full double-precision 
conversion of the program would increase the computing time 
by only 20% and would eliminate any problems of inaccuracy. 

4. SOME ADVANCES MADE FOR THE APPLICATION OF SEMI DIRECT 
METHODS 

We consider further now the use of direct elliptic solvers in 
semidirect methods and some further developments needed for 

1 
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applications in nonlinear fluid dynamics. Specifically con- 
sidered are the application to steady, inviscid transonic 
flow and the treatment of interior conditions. 

4.1 Transonic Flow . It is convenient to transform the 
perturbation velocities u, v {eqs, (2,4)), according to: 

u = (a/P)[l + u (x ,y) ] , v = av(x,y) , (4.1a,b) 

y — Y/P , X = x , P = (1 - M D l/t , a = ^ 3 /[ r ( T + (4.1c) 

Then eqs. (2.3), with f = -(l/2)u 2 , become 

u- + v- = u- - f- , 5 v- * 0 . (4.2a,b) 

x y x x y x 

For a nonlifting symmetrical thin airfoil, the chord can be 
located on one boundary (y = 0) , and the surface tangency 
condition can be approximated by 

v=av=F'(x) on y = 0 , (4.3) 

where y = <3 tf(x) is the equation of the airfoil surface 
between the leading and trailing edges, and* F(x)= 0 off the 
chord. Far from the airfoil, u * -1 and v -»■ 0 as x 2 + v 2 + 00 . 

If u < 0 everywhere, eqs. (4.2> are elliptic everywhere. 
But as M m increases from zero toward unity, both P and a 
become smaller, and then the boundary values of v become 
sufficiently large that u attains positive values in some 
regions, where eqs. (4.2) are then hyperbolic. Shock-wave 
discontinuities (in u) can exist where the flow decelerates 
from supersonic to subsonic velocity. Two related important 
factors are then the accuracy of a finite-difference repre- 
sentation and the iterative convergence of a scheme to solve 
the difference equations. 

For purely subsonic flow, (4.2) can be central differ- 
enced, and the semidirect method converges rapidly. But for 
transonic flow, the pure central differencing is not appro- 
priate in the embedded hyperbolic region, and the iteration 
with all central differences does not converge. Murman and 
Cole (see /76/) introduced type-dependent differencing in 
their line-relaxation method. In /76/, Murman refined the 
difference operators to make a fully conservative calculation, 
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which is important for capturing a shock wave in the proper 
location and for having the proper strength of the jump dis- 
continuity. Such a scheme uses upwind differencing in the 
hyperbolic regions to account for the domain of dependence; 
and transition operators for switching between the elliptic 
and hyperbolic operators provide the needed conservation 
property. To use type-dependent differencing in the semidi- 
rect method, note first that eq . (4.2a) is just -uu- + v- - 0, 
in which the first term is f-. Since the sign of u in f^ 
determines the type for eqs . (4.2), only f- needs type- 

dependent differencing. All other terms can be central 
differenced. With subscript C denoting central differences 
and subscript T denoting a type-dependent difference 
operator, a finite-difference representation of eqs. (4.2) is 

< 5 5>c + (5 y>c - < 3 ;>c - ' < 4 - 4al 

<5 y>C - <v 5 >c = 0 ■ (4 ' 4b> 

Although the use of type-dependent differencing ir _ eases 
the range of conditions (e.g., higher M^) for which the semi- 
direct iteration converges, it is not sufficient for flows 
with large supersonic zones. It is necessary /67,68/ to 
desymmetrize (with an "upwind bias") the iteration matrix 
representing the left side of eqs. (4.4). This can be done 
by adding a term, (-~a/Ax)u. to both the left and right 

1 1 r j 

sides of eq. (4,4a) for use in the iteration, where the 
central differences in (4.4a) are centered at point i',j' 
on Fig. 1. The effects on the stability and convergence of 
the iteration were explored in /67/, and appropriate values 
of 5 were determined in /68/. 

With simple central differences on a uniform mesh, all 
the terms other than (f-)m in eqs. (4.4) are second-order- 
accurate. The specification of (f -) T , however, determines 
(a) the formal order of accuracy of the solutions, (b) whether 
the desired conservation property is fulfilled, and (c) the 
capability of the difference equations to properly describe 
both the elliptic and hyperbolic portions of the flow field 
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and the transitions (such as through an abrupt "shock wave"). 

For the finite differences in eqs . (4.4), refer to Fig. 1. 

Equation (4.4b) is centered at the point i,j and the result 

is eq. (3.2b) with to. . = o. The shaded area on Fig. 1 is a 

mesh cell for the difference equation (4.4a), which represents 

the conservation equation f- + v- = 0. The central differ- 

x y 

ences in eq . (4.4a) are the same as in eq . (3.2a). Let 

(f x>T = (f \ V “ ! i-l ' ff 4. = “ (1/2) (u. -,) 2 . 

XT i ,3 1 1,3 1,3 1/3 { 4 . 5a , b) 

For central differencing, u, ., = u. for both f. and 

a 1,31/3 1/3 

f. For the general type-dependent differencing derived 

11,3 A 

in /68/, an expression for u. has been given in /69/ that 

i/3 

produces (f-) in eqs. (4.5) to either first- or second-order 

X X 

accuracy and automatically produces all four types of opera- 
tors (elliptic or hyperbolic and the transition operators, 
"parabolic" and "shock-point"). The result (dropping the 
primes from i‘ and j’ from here on) is 


u. . = (1 - a. . )u. .+ a. .[Xu. . + (1 - X)u. „.] , ( 

i/3 1/31/3 i/3 i -1 .3 '1-2,3 


4.6) 


where o . . is 0 or 1 depending on whether u. . < 0 or > 0, 

i/3 _ i/3 

that is, a. . = 0.5(1 + sgn u- .), where 
i/3 i/3 


+ “i, j’ ’ 


(4.7) 


and X is either 1 or (about) 2, corresponding to the formal 

order of accuracy desired. Note that u. . is effectively 

1,3 

the transformed velocity at the center of the shaded mesh cell 
but its sign (i.e., whether subsonic or supersonic) is used to 
determine how the flux f . . is represented at the downstream 

boundary of that mesh cell in (4.5a). "Artificial viscosity" 
in the second-order-accurate case is included by using 
X = ju + (1 - (i)/cAx, where M is 1 or 2 for first- or second- 
order accuracy and k is a positive number, 0(1) . The result 
from (4.5) is then 

(AX) (fj) T = - (1/2 ) j(l - o 1(j )5 i(j + 

+ u - MG..,, .]} 2 + d/2) j(l - 


+ cr. .[Xu. . + (1 — X)u. . ] | ^ 

1-1,3 1-2,3 i”3 , 3 ‘ 


(4.8) 
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i'‘jr use in eqs. (4.4) which become 

“i,j - 11 + + (A5/A ?> !a i,j - 

■ 5 i,j - (1 + - (A5 > (f 5»T ■ 

5 i,j+i ' 5 i,j - ( A 5/^Hv 1+lij - v i(j ) = 0 


(4.9a) 

(4.9b) 


In the semidirect iteration, extended direct Cauchy - 
Riemann solvers /74,71/ that can include the extra term (with 
coefficient a) solve the left side of eqs. (4.9) at all 
points in the field simultaneously at iteration n, with the 
right side evaluated at iteration n - 1. The four types of 
difference operators are implicitly included, and shock-wave 
discontinuities, represented by rapid changes over a few mesh 
points, are captured automatically. 


One can accelerate the iterative convergence by using 
shifting and scaling transformations /18/ and relaxation fac- 
tors (see /68/) and by a special technique for using the 
Aitken extrapolation formula /70/. 


Some results have been computed for a simple biconvex 
airfoil using an earlier inefficient version of this method. 
The outer conditions, at infinity, were replaced by prescrib- 
ing the "Prandtl-Glauert" solution on boundaries at 1/2-chord 
length fore and aft of the airfoil and at five chord lengths 
above the airfoil. Representative results for the thin- 
airfoil approximation to the surface pressure coefficient 
(C -2tu) over the chord are shown in Fig. 3. The mesh 

hr 

for Fig. 3(a) had 20 mesh points on the airfoil chord, and 
the results compare well with those from a line-relaxation 
program /77/, which used a refined variable mesh. All 
values of -C above C* are in a supersonic (hyperbolic) 

ir P 

region that is terminated downstream by a shock-wave dis- 
continuity. The shock is captured well by the finite- 
difference scheme. Fir .re 3(b) illustrates that the second- 
order-accurate results approximate the pressure distribu- 
tion well even on a very coarse mesh on which the first- 
order results are not accurate. 
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Fig. 3.- Pressure on a thin biconvex airfoil (M^ = 0.85, 

r = 0.10). 

The short computing time required is the most significant 
property of the semidirect method. The measured CPU time per 
iteration on a Control Data 7600 computer for the 39X32 mesh 
was 40 ms for an inefficiently coded program of an inefficient 
version of the algorithm. It is estimated for various 
reasons, discussed in /68/, that m e time can be reduced at 
least to 20 ms per iteration. The direct solver /74/ used 
14 ms on this mesh. The subcritical cases converged in 
3 iterations or less, and a slightly supercritical case 
(using an extrapolation) required only 6 iterations. The 
first-order-accurate case in Fig. 3(a) required 35 iterations. 
These methods are continually being improved, and more effi- 
cient computations are expected. 

4.2 In terior conditions . The above application to 
transonic flow was for a nonlifting biconvex airfoil- Because 
of symmetry, the airfoil chord could be located on one boundary 
of the computation field, and the airfoil surface condition 
could be transferred to y = 0. However, the real problem of 
interest in aerodynamics is the lifting airfoil (or wing) . 

Then there is generally no symmetry condition, and one needs 
to consider the airfoil in the interior of the flow field. 
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For treating interior conditions with direct elliptic 
solvers, the approaches can be grouped into two categories: 

(a) coordinate transformations, and (b) interior boundary- 
techniques without transformations. When conditions are to 
be applied on an interior surface such as on a two-dimensional 
airfoil, a transformation can put the interior surface on one 
boundary of the computation field of the transformed varjanles. 
Such transformations can be either analytical or numerical. 
These techniques are being vigorously pursued in other solu- 
tion methods, but have not been extensively employed with 
direct solvers, although Jameson (e.g,, /53/) has used an ana- 
lytical transofrmation when using a direct Poisson solver to 
accelerate SOR iterations. On the basis of work currently in 
progress, it is believed that general coordinate transforma- 
tions, including those done numerically, hold significant 
promise for use with direct solvers and semidirect methods. 

The other approach is to leave the geometry undistorted 
and deal with the interior conditions directly. Hockney /48/ 
adapted the classical capacity-matrix concept to numerical 
computation. That numerical technique, with all the subsequent 
refinements and generalizations that have been made for its 
applications, is referred to as the capacity-matrix technique 
(CMT) . In /66/ procedures were outlined for applying the 
technique with more general conditions than had been applied 
previously. The motivating physical problem was flow over a 
lifting airfoil. Specific procedures were given for applying 
conditions at points on an airfoil surface not coinciding with 
mesh points and for including the determination of the circu- 
lation (which produces the lift) about the airfoil to satisfy 
the Kutta condition at the trailing edge. In aerodynamics, 
that determination is normally done iteratively, but the 
"generalized" CMT provided the direct solution for the flow 
including the Kutta condition without iteration. One extra 
row and one extra column in the capacity matrix account for 
the Kutta condition and the circulation. 
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Although the examples used in /66/ were for Laplace's 
equation, the same technique has been used for calculation of 
nonlinear compressible flow over a lifting airfoil using a 
semidirect method for eqs • (2.2) in a preliminary version of 

/73/, with a condition = 0 on the airfoil surface. 

Recently various methods were studied / 1 2 / for treating 
conditions imposed on an interior slit. One method is the 
straightforward application of CMT with a Poisson operator. 
Extensions of this in a semidirect method can use combinations 
of CMT and shifting of matrix elements from one side of the 
iteration equation to the other. These techniques are useful 
for problems containing both a slit (where, e.g. , a Neumann 
condition is applied) and a cut behind the slit on which a 
discontinuity in the solution is enforced. In /72/ the CMT 
concept is also extended to use with first-order-elliptic 
equations, and various methods for applying discontinuous 
conditions on a slit are also explored. Such treatments can 
be useful in aerodynamic flows, in which either an airfoil 
(wing) can be approximated by a slit (sheet) or transforma- 
tions can be used that represent parts of airplanes (wings, 
fuselage, etc.) as sheets in a three-dimensional configuration. 

5 . CONCLUDING REMARKS 

The future prospects are encouraging for the use of direct 
elliptic solvers within semidirect methods for practical 
problems. On the basis of developments in progress it is 
anticipated that the full potential-flow equations (III, see 
section 2) will be treated successfully by semidirect methods 
for general inviscid transonic flows. Looking further into 
the future, based on some preliminary developments, it is 
believed that the most general fluid-dynamic equations (II 
and I) may eventually be solvable, for practical problems, 
using generalized forms of the semidirect method. 
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